*******************************************************************************
****						Table 5, panel A, Cols III-V - REGRESSION OF TOTAL TOLL ****
****		ON TIME DIFFERENTIALS: OTHER FUNCTIONAL FORM				*******
*******************************************************************************

use ".\data\clean\I10W_laneuse_dataset_15nov14_wcensus", clear

merge m:1 date hour using ".\data\clean\HV_ML_reliab.dta", keep(1 3) nogen 
gen reliabilityML=dist/p20_speedML-dist/p50_speedML
gen reliabilityHV=dist/p20_speedHV-dist/p50_speedHV
gen reliability_diff=reliabilityML-reliabilityHV
replace reliability_diff=0 if reliability_diff<0.01
keep if reliability_diff~=.
la var reliability_diff "Reliability"
*WTP calculation
gen TT_dif_hr=dist/MLspeed-dist/ELspeed
gen WTP2=charged_toll/TT_dif_hr
drop if ELspeed==.
drop if TT_dif_hr==.
drop if holiday==1
la var TT_dif_hr "Time in Hours"
drop if dow==0|dow==6
keep if acct_type=="PRIVATE"&occupancy~="HOV-3"
keep if hour>4 & hour<9
drop if TT_dif_hr<=0
forval p=1/5{
	g TT_dif_hr`p'=TT_dif_hr^`p'
	la var TT_dif_hr`p' "Time in Hours^`p'"
	}


estimates clear

*********************************************************
****	Col 1	Cubic Model
*********************************************************
eststo : reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3 reliability_diff  , cluster(rt_id)
// Root find to where WTP from TT saving==0
qui sum reliability_diff if e(sample)==1
scalar A0=_b[_cons]+`r(mean)'*_b[reliability_diff]
scalar A1=_b[TT_dif_hr]
scalar A2=_b[TT_dif_hr2]
scalar A3=_b[TT_dif_hr3]
forval p=1/3{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}

mata									
	a0 = st_numscalar("A0")
	a1 = st_numscalar("A1")
	a2 = st_numscalar("A2")
	a3 = st_numscalar("A3")
	s = polyroots((a0,a1,a2,a3))
	sReal = select(s, Im(s):==0)
	sReal = Re(sReal)
	sGeq0 = select(sReal, sReal[.,.]:>0)
	sLeq1 = select(sGeq0, sGeq0[.,.]:<1)
	sol = sLeq1
	st_numscalar("solution", sol)
	st_local("solution", strofreal(sol))
end
qui sum charged_toll if e(sample)==1
local toll1=round(`r(mean)',.01)
estadd scalar toll=`toll1'
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+ TT_dif_hr3*`mean3'
local wtp1=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp1'
local root1=round(60*`solution',.1)
estadd scalar root=`root1'

*********************************************************
****	Col 2	Quartic Model
*********************************************************

reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3 TT_dif_hr4 reliability_diff  , cluster(rt_id)
est sto a2

// Root find to where WTP from TT saving==0
forval p=1/4{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
qui sum charged_toll if e(sample)==1
local toll2=round(`r(mean)',.01)
estadd scalar toll=`toll2'	
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+ TT_dif_hr3*`mean3'+ TT_dif_hr4*`mean4'
local wtp2=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp2'


*********************************************************
****	Col 3	Quintic Model
*********************************************************

 
reg charged_toll TT_dif_hr TT_dif_hr2 TT_dif_hr3 TT_dif_hr4 TT_dif_hr5 reliability_diff  , cluster(rt_id)	
est sto a3
// Root find to where WTP from TT saving==0
qui sum reliability_diff if e(sample)==1
scalar A0=_b[_cons]+`r(mean)'*_b[reliability_diff]
scalar A1=_b[TT_dif_hr]
scalar A2=_b[TT_dif_hr2]
scalar A3=_b[TT_dif_hr3]
scalar A4=_b[TT_dif_hr4]
scalar A5=_b[TT_dif_hr5]
forval p=1/5{
	qui sum TT_dif_hr`p' if e(sample)==1
	local mean`p'=`r(mean)'
	}
qui sum charged_toll if e(sample)==1
local toll3=round(`r(mean)',.01)
estadd scalar toll=`toll3'	
lincom	TT_dif_hr*`mean1'+TT_dif_hr2*`mean2'+ TT_dif_hr3*`mean3'+ TT_dif_hr4*`mean4'+ TT_dif_hr5*`mean5'
local wtp3=round(`r(estimate)',0.01)
estadd scalar wtp=`wtp3'
mata									
	a0 = st_numscalar("A0")
	a1 = st_numscalar("A1")
	a2 = st_numscalar("A2")
	a3 = st_numscalar("A3")
	a4 = st_numscalar("A4")
	a5 = st_numscalar("A5")
	s = polyroots((a0,a1,a2,a3,a4,a5))
	sReal = select(s, Im(s):==0)
	sReal = Re(sReal)
	sGeq0 = select(sReal, sReal[.,.]:>0)
	sLeq1 = select(sGeq0, sGeq0[.,.]:<1)
	sol = sLeq1
	st_numscalar("solution", sol)
	st_local("solution", strofreal(sol))
end
local root3=round(60*`solution',.1)
estadd scalar root=`root3'
esttab *  using ".\results\appendix\tabs\Table5_pa_colIII_V.csv", replace  ///
	nonumbers nodepvars stats(N r2  aic bic ll wtp toll root, ///
		fmt(%9.0f %9.2f %9.0f %9.0f %9.0f %9.2f %9.2f %9.4f) ///
		labels(Observations R-Squared AIC BIC Log-Likelihood "WTP from Travel Time" "Total Tol" ///
		 "Time Savings in Min. for WTP=0")) ///
	cells(b(star fmt(%9.2f)) se(par)) star(* 0.10 ** 0.05 *** 0.01)  ///
	label order(_cons TT_dif_hr TT_dif_hr2 TT_dif_hr3 TT_dif_hr4 TT_dif_hr5  ///
		main:_cons b1:_cons b2:_cons b3: _cons reliability_diff)


